This document represents the Pytho codes associated with
Series Representations and Simulations of Isotropic Random Fields in the Euclidean Space
by Zhengwei Ma and Chunsheng Ma
# python imports
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import itertools
from numpy.random import normal, uniform, gamma
from scipy.special import jv
import plotly.graph_objects as go
import multiprocessing as mp
Here we set some initial values for simulation, as well as create some helper functions
r = 10
theta = 0.2
simulation_number = 1000 # number of times to simulate
m = 100 # number of points to simulate
r_list = np.linspace(0, 10, m)
theta_list = np.linspace(0, np.pi * 2, m)
R, T = np.meshgrid(r_list, theta_list) # R, T are matrices of size 100 by 100
X, Y = R*np.cos(T), R*np.sin(T) # express in polar coordinates
def generate_z(function):
z_list = []
with mp.Pool(mp.cpu_count()) as p:
z_list = p.starmap(function, itertools.product(range(m), range(m)))
Z = np.reshape(z_list, (m, m))
return Z
We want \begin{equation} y_{1} = -log(U_{0, 1}) \end{equation}
\begin{equation} y_{2} = U_{log(2), log(8)} \end{equation}\begin{equation} v_{1} = y_{1}^{0.5} * exp(\frac{-y_{2}}{2}) \end{equation}\begin{equation} w_{1} = U_{0, 1} \end{equation}\begin{equation} u_{1} = U_{0, 1} \end{equation}def example_1(i, j):
"""
This function takes in two positional arguments, i and j, which represents the corresponding
values in a m by m matrix (based on the value set in cell 4)
"""
r = R[i, j]
theta = T[i, j]
z_output = []
for i in range(simulation_number):
y1 = -np.log(uniform(low=0, high=1, size=1))
y2 = uniform(low=np.log(2), high=np.log(8), size=1)
v1 = (y1 ** 0.5) * np.exp(-y2/2)
w1 = uniform(low=0, high=1, size=1)
u1 = uniform(low=0, high=1, size=1)
Zi = jv(i, 2*r*np.sqrt(-v1 * np.log(w1))) * np.cos(i*theta+2*np.pi*u1)
z_output.append(Zi)
return np.sqrt(2) * np.sum(z_output)
Z = generate_z(example_1)
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_traces(contours_z=dict(
show=True, usecolormap=True,
highlightcolor="limegreen"))
fig.update_layout(title='', autosize=False,
width=1000, height=1000,
margin=dict(l=65, r=50, b=70, t=80), scene=dict(
xaxis=dict(showticklabels=False),
yaxis=dict(showticklabels=False),
zaxis=dict(showticklabels=False),
))
fig.show()
We want \begin{equation} y_{1} = \Gamma(0.5, 1) \end{equation}
\begin{equation} y_{2} = U_{2^{0.5}, 8^{0.5}} \end{equation}\begin{equation} v_{1} = \frac{y_{1}^{0.5}}{y_{2}} \end{equation}\begin{equation} w_{1} = U_{0, 1} \end{equation}\begin{equation} u_{1} = U_{0, 1} \end{equation}def example_2(i, j):
r = R[i, j]
theta = T[i, j]
z_output = []
for i in range(simulation_number):
y1 = gamma(shape=0.5, scale=1, size=1)
y2 = uniform(low=(2**0.5), high=(8**0.5), size=1)
v1 = (y1 ** 0.5) / y2
w1 = uniform(low=0, high=1, size=1)
u1 = uniform(low=0, high=1, size=1)
Zi = jv(i, 2*r*np.sqrt(-v1 * np.log(w1))) * np.cos(i*theta+2*np.pi*u1)
z_output.append(Zi)
return np.sqrt(2) * np.sum(z_output)
Z = generate_z(example_2)
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_traces(contours_z=dict(
show=True, usecolormap=True,
highlightcolor="limegreen"))
fig.update_layout(title='', autosize=False,
width=1000, height=1000,
margin=dict(l=50, r=40, b=80, t=75), scene=dict(
xaxis=dict(showticklabels=False),
yaxis=dict(showticklabels=False),
zaxis=dict(showticklabels=False),
))
fig.show()
We want \begin{equation} y_{1} = \Gamma(3, 1) \end{equation}
\begin{equation} y_{2} = U_{\frac{1}{64}, \frac{1}{4}} \end{equation}\begin{equation} v_{1} = y_{1}^{0.5} * y_{2}^{0.25} \end{equation}\begin{equation} w_{1} = U_{0, 1} \end{equation}\begin{equation} u_{1} = U_{0, 1} \end{equation}def example_3(i, j):
r = R[i, j]
theta = T[i, j]
z_output = []
for i in range(simulation_number):
y1 = gamma(shape=3, scale=1, size=1)
y2 = uniform(low=(1/64), high=(1/4), size=1)
v1 = (y1 ** 0.5) * (y2 ** 0.25)
w1 = uniform(low=0, high=1, size=1)
u1 = uniform(low=0, high=1, size=1)
Zi = jv(i, 2*r*np.sqrt(-v1 * np.log(w1))) * np.cos(i*theta+2*np.pi*u1)
z_output.append(Zi)
return np.sqrt(2) * np.sum(z_output)
Z = generate_z(example_3)
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_traces(contours_z=dict(
show=True, usecolormap=True,
highlightcolor="limegreen"))
fig.update_layout(title='', autosize=False,
width=1000, height=1000,
margin=dict(l=30, r=40, b=50, t=60), scene=dict(
xaxis=dict(showticklabels=False),
yaxis=dict(showticklabels=False),
zaxis=dict(showticklabels=False),
))
fig.show()
def besselnum(n):
results = []
i = 0
while i < n:
u = uniform(0, 1)
ug = uniform(0, 1)
g = 2 * ug/(1-ug)
if 3 * (0.5 / ((1+0.5*g) ** 2) * u) <= 2/g*(jv(g, 1) ** 2):
i+=1
results.append(g)
return results
def example_4(i, j):
r = R[i, j]
theta = T[i, j]
u = uniform(low=0, high=1, size=simulation_number)
v = besselnum(simulation_number)
z_list = []
for i in range(simulation_number):
Z = jv(i, r * v[i]) * np.cos(i*theta+2*np.pi*u[i])
z_list.append(Z)
return np.sqrt(2)*sum(z_list)
Z = generate_z(example_4)
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_traces(contours_z=dict(
show=True, usecolormap=True,
highlightcolor="limegreen"))
fig.update_layout(title='', autosize=False,
width=1000, height=1000,
margin=dict(l=50, r=50, b=55, t=80), scene=dict(
xaxis=dict(showticklabels=False),
yaxis=dict(showticklabels=False),
zaxis=dict(showticklabels=False),
))
fig.show()
For simulation of
\begin{equation} Z(x) = \sqrt(2) \sum_{n=0}^{\infty}J_{n}(rV_{n})cos(n\pi + 2\pi U_{n}) \end{equation}We being with simulating uniform random variables \begin{equation} U_{0}, U_{1}, U_{2} \sim [0,1]\end{equation}
rand_uniform = uniform(low=0, high=1, size=3) # simulate 3 random uniform
u0, u1, u2 = rand_uniform
The distribution function $ F(u) $ has to be chosen from say Table 1.
Consider the function on the second row of Table 1. When $d=2, \alpha=1$ we have \begin{equation} f(u) = \frac{1}{2} * u * exp(\frac{-u^2}{4}), u > 0 \end{equation} and \begin{equation} F(u) = 1 - exp(\frac{-u^2}{4}), u > 0 \end{equation}
Thus, we can use the inverse transformation method to generate random numbers from $F(u)$, where the inverse function of $F(u)$ is \begin{equation} F^{-1}(u) = 2 * \sqrt{-\log(1-u)} \end{equation}
For a uniform random variable $U$ on (0, 1), $U$ and $1-U$ have the same distribution. Thus, a code to generate one $F(u)$ random number is
\begin{equation}v = 2 * \sqrt{\log{U_{0}}} \end{equation}rand_v_uniform = uniform(low=0, high=1, size=3)
v0, v1, v2 = 2 * np.sqrt(-np.log(rand_v_uniform))
Demonstrate Z(x) in formula 3.7 with an example, with n=2
np.sqrt(2) * (
jv(0, r*v0) * np.cos(0*theta+2*np.pi*u0) + \
jv(1, r*v1) * np.cos(1*theta+2*np.pi*u1) + \
jv(2, r*v2) * np.cos(2*theta+2*np.pi*u2)
)
-0.05915042467644759
We show that the following formula is equivalent
def func_v(x):
return 2 * np.sqrt(-1 * np.log(x))
np.sqrt(2) * sum(
[
jv(i, r*func_v(rand_v_uniform[i])) * \
np.cos(i*theta+2*np.pi*rand_uniform[i])
for i in range(0, 3, 1)
]
)
-0.05915042467644759
def generate_function(n, r, theta):
n_random_uniform = uniform(size=(n+1))
n_random_v_uniform = uniform(size=(n+1))
return np.sqrt(2) * sum(
[
jv(i, r*func_v(n_random_v_uniform[i])) * \
np.cos(i*theta+2*np.pi*n_random_uniform[i])
for i in range(0, n+1, 1)
]
)
# now calculate z for each r, theta
z_list = []
for theta in theta_list:
for r in r_list:
z_list.append(generate_function(200, r, theta))
## the output here is m_theta by m_r
Z = np.reshape(z_list, (m, m))
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_traces(contours_z=dict(
show=True, usecolormap=True,
highlightcolor="limegreen",
project_z=True))
fig.update_layout(title='Isotrophic Random Fields', autosize=False,
width=1000, height=1000,
margin=dict(l=60, r=50, b=60, t=90), scene=dict(
xaxis=dict(showticklabels=False),
yaxis=dict(showticklabels=False),
zaxis=dict(showticklabels=False),
))
fig.show()